RRSVD background adjustment

Matrix Decomposition Method, Full Data

Vo and Park (2016) used a matrix decomposition approach to estimate background, foreground and noise in EM images. The math behind their method was worked out by Zhang et al. (2013) That is, if \(Y\) is a \(m\times n\) matrix where \(y_{ij}\) is the intensity value at the \(x\) and \(y\) locations. Then \(Y\) can be decomposed into \(B\), \(F\) and \(E\) corresponding to background, foreground and error, respectively, through \[Y = B + F + E.\]

\(B\) is estimated first using a regularized low-rank approximation \[\hat B_\lambda=\sum_{k=1}^p\hat{\mathbf{u}}_\lambda\hat{\mathbf{v}}_\lambda^\top\] where \(\hat{\mathbf{u}}_\lambda\) is the \(m\times 1\) unnormalized left singular vector estimated using regularized, robust SVD with penalty parameter \(\lambda\) and \(\hat{\mathbf{v}}_\lambda\) is the \(n\times 1\) right singular vector with unit norm estimated using RRSVD with penalty \(\lambda\). The rank of the approximation and the penalty parameter \(\lambda\) needs to be chosen on a case by case basis. A few example estimates of \(B\) with different \(\lambda\) values is below. Each of the estimates is an \(8^{th}\) degree approximation. Similar to regularized regression, large values of \(\lambda\) smooth out the background estimate while small values result in a very s.

#Create 'Bhat_list' object, takes a while so only run once and load the result in the next chunck
#Robust SVD of ex_frame in order to estimate background matrix (section 2.1), rank p
p <- 1
lams <- c(0,1,10,20,50)
Bhat_list <- vector("list",length(lams))
s_s <- matrix(0,length(lams),p)
for(j in 1:length(lams)){
  Bhat <- matrix(0,256,256)
  y <- ex_frame
  for(i in 1:p){
    y <- ex_frame - Bhat
    rsvdi <-  RobRSVD(y, irobust=TRUE, niter=100,uspar = lams[j], vspar = lams[j])
    Bhat <- Bhat + rsvdi$u%*%t(rsvdi$v)*rsvdi$s
    s_s[j,i] <- rsvdi$s
  }
  Bhat_list[[j]] <- Bhat
}
#plot(s_s[2,-c(1:2)])
#save(Bhat_list,file="/Users/stan070/Documents/AIM Project (Lisa)/streaming-tcs/EM_Use_Case/Write_Ups/Bhat.Rdata")

load("~/Documents/AIM Project (Lisa)/streaming-tcs/EM_Use_Case/Write_Ups/Bhat.Rdata")
Vo_method_df <- data.frame(x=m1_dat$x,y=m1_dat$y, Y = matrix(t(ex_frame),ncol=1))
Vo_method_df$Bhat <- matrix(t(Bhat_list[[1]]),ncol=1)
plot1 <- qplot(x,y,colour=Bhat,data=Vo_method_df)+theme_em()+ggtitle(expression(lambda==0))
Vo_method_df$Bhat <- matrix(t(Bhat_list[[2]]),ncol=1)
plot2 <- qplot(x,y,colour=Bhat,data=Vo_method_df)+theme_em()+ggtitle(expression(lambda==1))
Vo_method_df$Bhat <- matrix(t(Bhat_list[[3]]),ncol=1)
plot3 <- qplot(x,y,colour=Bhat,data=Vo_method_df)+theme_em()+ggtitle(expression(lambda==10))
Vo_method_df$Bhat <- matrix(t(Bhat_list[[4]]),ncol=1)
plot4 <- qplot(x,y,colour=Bhat,data=Vo_method_df)+theme_em()+ggtitle(expression(lambda==20))
Vo_method_df$Bhat <- matrix(t(Bhat_list[[5]]),ncol=1)
plot5 <- qplot(x,y,colour=Bhat,data=Vo_method_df)+theme_em()+ggtitle(expression(lambda==50))
grid.arrange(plot1,plot2,plot3,plot4,plot5,nrow=1)

Vo_method_df$Bhat <- matrix(t(Bhat_list[[2]]),ncol=1)
rm(Bhat_list)

The decomposition of the original data (\(Y\), left plot) into background (\(B\), middle plot) and leftover (\(F+E\), right plot) using \(\lambda=1\) is below. The right plot does seem to have the left-to-right gradient removed, but it lightens the node up a bit too much.

Vo_method_df$Ohat <- Vo_method_df$Y - Vo_method_df$Bhat
truth <- qplot(x,y,colour=Y,data=Vo_method_df)+theme_em()+ggtitle("Y (original)")
bground <- qplot(x,y,colour=Bhat,data=Vo_method_df)+theme_em()+ggtitle("B (background)")
leftover <- qplot(x,y,colour=Ohat,data=Vo_method_df)+theme_em()+ggtitle("F+E (leftover)")
grid.arrange(truth,bground,leftover,nrow=1)

p <- 16
lambda <- 1
imgs <- vector("list",p)
all_plots <- vector("list",p)
approx <- matrix(0,256,256)
s_s <- rep(NA,p)
for(i in 1:p){
  y <- ex_frame - approx
  rsvdi <-  RobRSVD(y, irobust=TRUE, niter=500,uspar = lambda, vspar = lambda)
  imgs[[i]] <-  rsvdi$u%*%t(rsvdi$v)*rsvdi$s
  approx <- approx + imgs[[i]]
  s_s[i] <- rsvdi$s
  one_d_df <- data.frame(x=m1_dat$x,y=m1_dat$y, Int=matrix(t(imgs[[i]]),ncol=1))
  all_plots[[i]] <- qplot(x,y,colour=Int,data=one_d_df)+theme_em()+ggtitle(paste("Dimension",i))
}
#plot(s_s)
grid.arrange(grobs=all_plots,nrow=4,ncol=4)

#sum(s_s)/(256^2)

They then use the robust, regularized SVD approach again to decompose what’s left over (\(O\)) into foreground and residual, but I haven’t gotten that to give consistent results. Instead, I apply Maggie’s EM approach with no background to cluster the “left-over” pixels into growth and non-growth.

#Apply EM without background adjustment
mean_res <- EM_alg_mu_cpp(z=Vo_method_df$Ohat,p = rep(1/5,5), mu = c(-50,0,25,50,100), sigma = rep(4,5),eps = 1e-5,maxit = 1000)
Vo_method_df$Mean_class <- as.factor(apply(mean_res$wmat,1,which.max))
g_nong <- qplot(x,y,colour=Mean_class,data=Vo_method_df,size=I(.1))+theme_em()+scale_colour_manual(values=c(1,1,2,2,2))+ggtitle("Growth")
grid.arrange(truth,g_nong,nrow=1)

#hist(Vo_method_df$Ohat,breaks=1000)

Matrix decomposision method, sub-image

Now to do the same thing but in the presence of missing data. The SVD approach does not work for missing data so the missing data must be imputed before the SVD approach can be applied. Zhang et al. (2013) suggest the missing data be first replaced by the column or row mean, then estimate \(u\) and \(v\), replace missing data with \(u_iv_j\) and repeat until convergence, i.e., the estimates of \(u\) and \(v\) don’t change much. They say that starting the algorithm with the row or column mean makes no difference, but in our case the node screws everything us. Using the column means results in a band of dar columns in the middle and two lighter columns on the side. The row means come are very dark for the bottom half of the image and much lighter on the top. Therefore I impute the missing cells with the mean of the column and row means, which appears to work okay. The spline imputation method Maggie and Sarah have been using may provide better results but I haven’t tried yet.

subset <- m1_dat[,c("x","y","T2621")]
subset$T2621[sample(nrow(m1_dat),size=(.9*nrow(m1_dat)))] <- NA
ex_frame_samp <- matrix(subset$T2621,nrow=256,ncol=256,byrow = TRUE)
Bhat_missing <- lr_approx_missing(Y = ex_frame_samp,p = 1,lambda = 1)
Vo_method_df$Bhat_missing <- matrix(t(Bhat_missing$approx),ncol=1)
Vo_method_df$Ohat_missing <- Vo_method_df$Y - Vo_method_df$Bhat_missing
truth <- qplot(x,y,colour=Y,data=Vo_method_df)+theme_em()+ggtitle("Y (original)")
bground <- qplot(x,y,colour=Bhat_missing,data=Vo_method_df,size=I(.2))+theme_em()+ggtitle("B (background)")
leftover <- qplot(x,y,colour=Ohat_missing,data=Vo_method_df,size=I(.2))+theme_em()+ggtitle("F+E (leftover)")
grid.arrange(truth,bground,leftover,nrow=1)

Apply Maggie’s EM algorithm to the background corrected image. There results are for five groups. The red dots are groups 4 and 5 (which are considered growth) while the green dots are groups 1-3.

no_bg <- Vo_method_df$Ohat_missing
no_bg <- no_bg[-which(is.na(no_bg))]
#hist(no_bg,breaks=500)
mean_res_missing <- EM_alg_mu_cpp(z=no_bg,p = rep(1/5,5), mu = c(-50,-10,0,50,100), sigma = rep(4,5),eps = 1e-5,maxit = 1000)
Vo_method_df$Mean_class_missing <- Vo_method_df$Ohat_missing
Vo_method_df$Mean_class_missing[which(!is.na(Vo_method_df$Mean_class_missing ))] <-   apply(mean_res_missing$wmat,1,which.max)
non_missing_df <- subset(Vo_method_df,!is.na(Vo_method_df$Mean_class_missing))
non_missing_df1 <- subset(non_missing_df,Mean_class_missing>3)
non_missing_df2 <- subset(non_missing_df,Mean_class_missing<=3)
g_nong_missing <- ggplot(data=Vo_method_df)+geom_point(aes(x,y,colour=Y))+theme_em()+ggtitle("Estimated Growth") + geom_point(aes(x,y),colour=2,data=non_missing_df1)
g_nong_all <- ggplot(data=Vo_method_df)+geom_point(aes(x,y,colour=Y))+theme_em()+ggtitle("Estimated Growth") + geom_point(aes(x,y),colour=2,data=non_missing_df1)+
  geom_point(aes(x,y),colour='green',data=non_missing_df2)
grid.arrange(truth,g_nong_missing,g_nong_all,nrow=1)

Matrix Decomposition to Identify Growth Beginning

The singular values are a measure of how much of the signal in the image can be explained by the corresponding rank 1 approximation to the image. One would think that rank one matrices will do a better job of representing images without growth, i.e., only have a node and background, compare to images where there is a lot of growth going on. Test this by plotting the first singular value for the frames over time.

How about some of the later rank one matrices? Plot the singular values for the first five rank one approximations for five frames that exhibit no growth, initial stages of growth, max growth and dissipation.

try <- c("T199","T2162","T2287","T2621","T3666")
subset <- m1_dat[,c("x","y",try)]
s_mat <- matrix(0,length(try),5)
for(i in 3:ncol(subset)){
  subset[sample(nrow(m1_dat),size=(.9*nrow(m1_dat))),i] <- NA
  ex_frame_samp <- matrix(subset[,i],nrow=256,ncol=256,byrow = TRUE)
  Bhat_missing <- lr_approx_missing(Y = ex_frame_samp,p = 5,lambda = 1)
  s_mat[,(i-2)] <- Bhat_missing$ss
}
s_mat <- data.frame(s_mat)
s_mat_melt <- reshape2::melt(s_mat)
s_mat_melt$Index <- rep(1:5,5)

frame1 <- qplot(x,y,colour=T199,data=m1_dat)+theme_em()+ggtitle("X1")
frame2 <- qplot(x,y,colour=T2162,data=m1_dat)+theme_em()+ggtitle("X2")
frame3 <- qplot(x,y,colour=T2287,data=m1_dat)+theme_em()+ggtitle("X3")
frame4 <- qplot(x,y,colour=T2621,data=m1_dat)+theme_em()+ggtitle("X4")
frame5 <- qplot(x,y,colour=T3666,data=m1_dat)+theme_em()+ggtitle("X5")


full_pic <- qplot(Index,value,data=s_mat_melt,group=variable,colour=variable,geom='line')+geom_point()+theme(legend.position='none')+theme_bw()+
  ylab("Singular Value")+xlab("Dimension")
red_pic <- qplot(Index,value,data=s_mat_melt,group=variable,colour=variable,geom='line')+geom_point()+ylim(c(0,2000))+theme_bw()+xlim(2,5)+
  ylab("Singular Value")+xlab("Dimension")

lay <- rbind(c(1,2,3,4,5),
             c(6,6,7,7,7))

grid.arrange(frame1,frame2,frame3,frame4,frame5,full_pic,red_pic,layout_matrix=lay)